function [e,gamma] = computeGammaOpening(element,d,e0)

%code 1
% A = element.A;
% H = element.H;
% Dnn = element.Dnn;
% Dss = element.Dss;
% Dhat = [Dnn 0;0 Dss];
% ef = element.ef;
% e = e0;
% 
% % A*d
% % A*H
% % ef
% % Dnn
% % element.ttilde
% 
% etilde = computeEquivalentSeparation(element,e);
% gammas = 1/computeDamageFunction(element,etilde);
% if etilde < ef
%     Fs = gammas;
%     dFdetilde = computeDamageFunctionDerivative(element,etilde);
%     f = dFdetilde*Dhat*e/Dnn/etilde;
%     fs = -f*gammas^2;
%     lhs = gammas*Dhat + A*H + Dhat*e*fs';
%     rhs = A*(d-H*e) - Fs*Dhat*e;
% else
%     lhs = A*H;
%     rhs = A*(d-H*e);
% end
% 
% conv_criteria = 0;
% while conv_criteria == 0
%     de = lhs\rhs;
%     e = e+de;
%     
%     etilde = computeEquivalentSeparation(element,e);
%     gammas = 1/computeDamageFunction(element,etilde);
%     if etilde < ef
%         Fs = gammas;
%         dFdetilde = computeDamageFunctionDerivative(element,etilde);
%         f = dFdetilde*Dhat*e/Dnn/etilde;
%         fs = -f*gammas^2;
%         lhs = gammas*Dhat + A*H + Dhat*e*fs';
%         rhs = A*(d-H*e) - Fs*Dhat*e;
%     else
%         lhs = A*H;
%         rhs = A*(d-H*e);
%     end
%     
%     if norm(rhs) < 1e-9
%         conv_criteria = 1;
%     end
% end
% gamma = 1/gammas;

%code 2
e = e0;
etilde = computeEquivalentSeparation(element,e);
gammas = 1/computeDamageFunction(element,etilde);
A = element.A;
H = element.H;
Dnn = element.Dnn;
Dss = element.Dss;
Dhat = [Dnn 0;0 Dss];
conv_criteria = 0;
while conv_criteria == 0
    e = (gammas*Dhat + A*H)\(A*d);
    etilde = computeEquivalentSeparation(element,e);
    gammas = 1/computeDamageFunction(element,etilde);
    rhs = A*d - (gammas*Dhat + A*H)*e;
    if norm(rhs) < 1e-14
        conv_criteria = 1;
    end
end
gamma = 1/gammas;

end
